APPROXIMATION OF THE DETERMINANT OF LARGE SPARSE 
SYMMETRIC POSITIVE DEFINITE MATRICES 



ARNOLD REUSKEN* 

Abstract. This paper is concerned with the problem of approximating det{A)^^" for a large 
sparse symmetric positive definite matrix A of order n. It is shown that an efficient solution of 
this problem is obtained by using a sparse approximate inverse of A. The method is explained 
and theoretical properties are discussed. A posteriori error estimation techniques are presented. 
" Furthermore, results of numerical experiments are given which illustrate the performance of this new 

^rr^ ' method. 
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, 1. Introduction. Throughout this paper, A denotes a real symmetric positive 

' definite matrix of order n with eigenvalues 

^ ! 

< Ai < A2 < ... < A„. 

^ , In a number of apphcations, for example in lattice Quantum Chromodynamics |l^ , 

' certain functions of the determinant of A, such as det(^)^/" or ln(det(A)) are of 

, interest. It is well-known (cf. also ^) that for large n the function A det{A) has 

' poor scaling properties and can be very ill-conditioned for certain matrices A. In this 

, paper we consider the function 

o ■ 

O ; d: A~>det{A)i . (1.1) 

^ ' A few basic properties of this function are discussed in In this paper we present 

' . a new method for approximating d{A) for large sparse matrices A. The method is 

^ based on replacing A by a matrix which is in a certain sense close to A~^ and for 

. which the determinant can be computed with low computational costs. One popular 

^ ' method for approximating A is based on the construction of an incomplete Cholesky 

• i-H . factorization. This incomplete factorization is often used as a preconditioner when 

' solving linear systems with matrix A. In this paper we use another preconditioning 

^ . technique, namely that of sparse approximate inverses (cf. [|l], ^ In Re- 



mark 3.1c we comment on the advantages of the use of sparse approximate inverse 
preconditoning for approximating d{A). Let A = LL'^ be the Cholesky decomposi- 
tion of A. Then using techniques known from the literature a sparse approximate 
inverse Ge of L, i.e. a lower triangular matrix Ge which has a prescribed sparsity 
structure E and which is an approximation of i^^, can be constructed. We then 
use det(G_B)~^/" — n"=i(^B)ii^^" approximation for d(A). In §|| we explain 

the construction of G^; and discuss theoretical properties of this sparse approximate 
inverse. For example, such a sparse approximate inverse can be shown to exist for 
any symmetric positive definite A and has an interesting optimality property re- 
lated to d{A). As a direct consequence of this optimality property one obtains that 
d(A) < det(G£;)"2/" holds and that the approximation of d{A) by det(G£;)"2/" be- 
comes better if a larger sparsity pattern E is used. 
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In §^ we consider the topic of error estimation. In the paper j^] bounds for the deter- 
minant of symmetric positive definite matrices are derived. These bounds, in which 
the Frobenius norm and an estimate of the extreme eigenvalues of the matrix involved 
are used, often yield rather poor estimates of the determinant (cf. experiments in (P). 



In §4.1 we apply this technique to the preconditioned matrix Ge-AG"^ and thus obtain 



reliable but rather pessimistic error bounds. It turns out that this error estimation 



technique is rather costly. In §4.2 we introduce a simple and cheap Monte Carlo tech- 
nique for error estimation. In §^ we apply the new method to a few examples of large 
sparse symmetric positive definite matrices. 

2. Preliminaries. In this section we discuss a few elementary properties of the 
function d. We give a comparision between the conditioning of the function d and 
of the fuction A d{A)"- — det{A). We use the notation || ■ ||2 for the Euchdean 
norm and k{A) — A„/Ai denotes the spectral condition number of A. The trace of 
the matrix A is denoted by tr{A). 

Lemma 2.1. Let A and 6 A be symmetric positive definite matrices of order n. 
The following inequalities hold: 

Ai < d{A) < A„ , (2.1a) 

d{A) < -tT{A) , (2.1b) 
n 

diA + SA)-d{A) \\SAh 
° < d(A) - "^^^ \\Ah ■ ^ ' 



Proof. The result in (2.1a) follows from 

n 

Ai < (H^*)" 



The result in ( 2.1b ) follows from the inequality between the geometric and arithmetic 
mean: 



n n 

d{A) = (nX,)^ <_Va, = -tr(A) . 

J. A 77 ^ 97 



n 

1 = 1 



From the Courant-Fischer characterization of eigenvalues it follows that 

X^{A + 6 A) > A, (A) + Ai(M) > A, (A) 
for aU i. Hence d{A + SA) > d{A) holds. Now note that 

^iA±M_^.(det(/ + A-M))^-l 

= (j[il + \M~'5A))^ I 

^ (n(i+ii^"'ii2ii'^^ii2)) -1 

^\\A-'MSAh = '^(^)f^- 

ll^lb 
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Thus the resuh in ( 2.1c ) is proved. □ 

The resuh in ( 2.1c ) shows that the function d{A) is weU-conditioned for matrices 
A which have a not too large condition number k(A). 

We now briefly discuss the difference in conditioning between the functions A — > 
d{A) and A —>■ det{A). For any symmetric positive definite matrix B of order n we 
have 

diA + tB)-diA) ^ <MtriA-^B) . 
^ ^ t-*o t n ^ ' 

From the Courant-Fischer eigenvalue characterization we obtain Ai(A^^i?) < Ai(A^^)||_B| 
for all i. Hence 

„„.xM, \d'{A)B\ d{A) trM-ifi) d{A) ,^ 

M M) 2 max ' , = max \ „ — < -^trM^^) , 

" ^ '^"^ BisSPD ||B||2 n BisSPD ||B||2 " n ^ 

with equality for B = I. Thus for the condition number of the function d we have 
\\AUd'{A)h 1 



d{A) 



-\\A\\MA-^) < <A) . (2.2) 



Note that for the diagonal matrix A = diag(Aii) with An = 1, An = a G (0, 1) for 
i > 1, in the inequality in (^.2[) one obtains equality for n ^ oo. For this A and with 



5A = el, e > 0, for n — > cxd we have equality in the second inequality in (2.1c), too. 
For d{A) — det(A) — (i(A)" the condition number is given by 

d{A) d{AY II IM V ; , V ; 

i.e. n times larger than the condition number in ( p.2| ). The condition numbers for 
d and d give an indication of t he se nsitivity if the perturbation ||i5j4||2 is sufficiently 



small. Note that the bound in ( 2.1c ) is valid for arbitrary symmetric positive definite 
perturbations 5 A. The bound shows that even for larger perturbations the function 
d{A) is well-conditioned at A if k.{A) is not too large. For the function d{A) the 
effect of relatively large perturbations can be much worse than for the asymptotic 



case {5A 0), which is characterized by the condition number in (2.3). Consider, for 
example, for < e < ^ a perturbation 5A = eA, i.e. ||(5A||2/|| A||2 = e. Then 

^(A±ai^^(l + e)"-l>e>-l, 

d{A) y ) - 

which is very large if, for example, e = 10~^, n = 10^. 

The results in this section show that the numerical approximation of the function 
d{A) is a much easier task than the numerical approximation of A — > det{A). 

3. Sparse approximate inverse. In this section we explain and analyze the 
construction of a sparse approximate inverse of the matrix A. Let A = LL"'" be the 
Cholesky factorization of A, i.e. L is lower triangular and L~^AL~^ — I. Note that 
d{A) — d{LY ~ nr=i ^i/"- We will construct a sparse lower triangular approximation 
G of and approximate d{A) by d{G)^'^ — YVi=i^ii^^'^ ■ The construction of a 
sparse approximate inverse that we use in this paper was introduced in ^ ^] and 
can also be found in Q. Some of the results derived in this section are presented in 
0, too. 
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We first introduce some notation. Let E C | 1 < < n} be a given 

sparsity pattern. By we denote the number of elements in E. Let Se be tlie set 
of n X n matrices for which aU entries are set to zero if the corresponding index is not 
in E: 

Se = {M e R"""" I My = if (ij) i E} . 

YoY \ < i < n XnX Ei = E r\ {(i, j) | 1 < J < n}. If rn := =f^Ei > we use the 
representation 



< jn^ < n . (3.1) 
(3.2) 



Et = {(i,jl),(«,j2),--- ,(«,jni)}, 1 < jl <j2 < 

For rii > we define the projection 
Note that the matrix 

P,AP^ : R"' ^ R"' 

is symmetric positive definite. Typical choices of the sparsity pattern E (cf. S|^) are 
such that rii is a very small number compared to n (e.g. < 20). In such a case the 
projected matrix PiAP^ has a small dimension. 

To facilitate the analysis below, we first discuss the construction of an approx- 
imate sparse inverse Me G 5'_b in a general framework. For AIe 6 Se we use the 
representation 



nii e 



Note that if = then 
For given A,B gR" 
lowing problem: 



' = (0,0,... ,0). 

with A symmetric positive definite we consider the fol- 



determine Me £ Se such that {MEA)ij — Bij for all G E 



(3.3) 



In ( |3.3| ) we have #i? equations to determine #i? entries in Me- We first give two basic 
lemmas which will play an important role in the analysis of the sparse approximate 
inverse that will be defined in (3.9). 

Lemma 3.1. The problem ^3.^ ) has a unique solution Me G Se- IJ ni > then 
the ith row of Me is given by mf with 



TO, = Pi {P,AP^)-^PA , 

where bf is the ith row of B. 

Proof The equations in (|3.3[) can be represented as 



(3.4) 



{mf A)ji^ — {bj)j^ for all i with > and all fc = 1, 2, . . . , , 

where mf is the ith row of Me- Consider an i with > 0. Note that Me G Se, hence 
Pj^ Pimi = mi- For the unknown entries in mi we obtain the system of equations 



{AP,Pf^7n, 



1,2, 
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which is equivalent to 



The matrix PiAP^ is symmetric positive definite and thus rm must satisfy 

P,m, = {P,APl)-^P,b, . 



Using Pj" Piirii = we obtain the result in (3.4). The construction in this proof 
shows that the solution is unique. □ 

Below we use the Frobenius norm, denoted by || • \\f'- 



\B\\l=Y.Bl=t,{BB^) , Be 



(3.5) 



Lemma 3.2. Let A = LL^ be the Cholesky factorization of A and let Me G Se 
be the unique solution of ^3.^ ). Then Me is the unique minimizer of the functional 

M ^\\{B - MA)L-'^fp = tr((S - MA)A-^{B - MAf), M e Se- 

(3.6) 



Proof. Let be the ith basis vector in R". Take M € Se- The zth rows of M 
and B are denoted by mf and bf , respectively. Now note 

n 

tr((B - MA)A-^{B - MAf ) ^ ^ eJiBA^'^B'^ - MB^ - BM^ + MAM^)ei 

i=l 

n 

= tr{BA-'^B^) + ^(-2mf fe; + mjAmi) . (3. 7) 



The minimum of the functional ( 3.6 ) is obtained if in ( p.7| ) we minimize the functionals 

m, -2mfb, + mf Ami , rui € 7^(i^^) (3.8) 

f or a ll i with > 0. If we write m.^ = P^rhi , nii G M"' , then for > the functional 
(3.8) can be rewritten as 

rhi -2mfPA + mJPiAPj'mi , rhi G M"'. 

The unique minimum of this functional is obtained for rhi — {PiAP^)^^ Pibi, i.e. 
m^ = Pf{P,APf)-'^Pibi for aU i with > 0. Using Lemma ^ it follows that Me 
is the unique minimizer of the functional (3.6). □ 



Sparse approximate inverse. We now introduce the sparse approximate inverse 
that will be used as an approximation for L^^ . For this we choose a lower triangular 
pattern C {(j,j) | ^ < j < i < n} and we assume that {i,i) G for all i. The 
sparse approximate inverse is constructed in two steps: 

1. Gei e Sei such that {GEiA)ij = 5ij for all € E^ , (3.9a) 



2. Gb, (diag(G£0)"^G£, 



(3.9b) 
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The construction of G^i in (3_^) was first introduced in A theoretical background 
for this factorized sparse inverse is given in [pl| . The approximate inverse G^i in 
( 3.9a ) is of the form ( |3.3| ) with B = I. From Lemma 3T it follows that in ( 3.9a ) there 
is a unique solution G^i . Note that because is lower triangular and (i,i) G we 
have Hi = 4j=E^ > for all i and j„. = i in (3.1). Hence it follows from Lemma 3.1 
that the ith row of G^i , denoted by gf , is given by 



Pf{P,AP^yH,, with e, = (0, 



,0,1) 



(3.10) 



The zth entry of gi, i.e. 



T 



is given by ej {PtAP^ ) Si, which is strictly positive 



because PiAP^ is symmetric positi ve def inite. Hence diag(G'£;!) contains only strictly 
positive entries and the secon d ste p ( 3.9b ) is well-defined. Define cji = Pigi- The sparse 
approximate inverse G^ji in ( 3.9a ) can be computed by solving the (low dimensional) 
symmetric positive definite systems 



P,AP^g,^{0,... 



1,2, 



(3.11) 



We now derive some interesting properties of the sparse approximate inverse as in 
(3.9). We start with a minimization property of G^r. 

Theorem 3.3. Let A = LL?" he the Cholesky factorization of A and D := 
diag(L), L := LD. G^i as in (3.9a) is the unique minimizer of the functional 



G ||(/ - GL)D-^fp = tr((/ - GL)D-^{I - GLf), GeSe^. 



(3.12) 



Proof. The construction of G^' ii^ (3.9a) is as in ( |3.3| ) with E — E^ , B = I. Hence 
Lemma 3.2 is applicable with B = I. It follows that G^i is the unique minimizer of 



\\{I-GA)L- 



-T\\2 

\\f 



G E S j^i 



(3.13) 



Decomposed as L ^ = D ^ -l-i? with i? strictly upper triangular. Then _D ^—GL 
and R are lower and strictly upper triangular, respectively, and we obtain: 



||(/ - GA)L-^\\l =\\{I- GLL^)L-^\\l = \\D-' 



= \\D- 



GL\\% 



\R\\% = - GL)D 



R-GLfp 



-1 l|2 
F 



\R\\ 



Hence the minimizers in ( 3.13 ) and ( 3.12 ) are the same. 



□ 



Remark 3.4. From the result in Theorem 3.3 we see that in a scaled Frobenius 
norm (scaling with D^^) Gpi is the optimal approximation of L^^ in the set Spi, in 
the sense that Gpi L is closest to the identity. A seemingly more natural minimization 
problem is 



min 11/ 



GLW 



(3.14) 



i.e. we directly approximate L ^ (ins tead of L ^) and do not use the scaling with 
. The minimization problem (3.14) is of the form as in Lemma 3.2 with B = , 
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E = EK Hence the unique minimizer in ( 3.14 ), denoted by Ge'i must satisfy (3.3) 
with B ^ L"^: 



(G£;i^),y = Lj^ for aU (ij) e E^ 



(3.15) 



Because E^ contains only indices with i > j and Lji — for i > j it foUows that 
G^i & Se' rnust satisfy 



if*^J 
La if i = j 



for all (ij) G E^ 



(3.16) 



This is similar to the system of equations in (3.9a), which characterizes G^i ■ However, 
in ( 3.16D one needs the values La, which in general ar e not available. Hence opposite 
to the minimization problem related to the functional (3.12) the minimization problem 
( 3.14 ) is in general not solvable with acceptable computational costs. □ 

The following lemma will be used in the proof of Theorem 3.7. 

Lemma 3.5. Let Ge' be as in (3.9a). Decompose Ge' as Ge' = D{I — L), with 
D diagonal and L strictly lower triangular. Define E[_ E^ \ {{i.,i) | 1 < * < n}. 
Then L is the unique minimizer of the functional 



L tr((/ - L)A{I - L^)) , LeSei , 
and also of the functional 

L det[diag((/ - L)A{I ~ L'^))] , Le Se 

Furthermore, for D we have 

D = [diag((/ - L)A{I - L^))]-i . 



(3.17) 



(3.18) 



(3.19) 



Proof. From the construction in (3.9a) it follows that 
((/ - L)A),j = for all {i,j) £ E[ 



i.e., L e Se' is such that {LA)ij — Aij for all (i, j) G 5^;! . This is of the form (3.3) 
with B = A^ E — E^_. From Lemma 3.2 we obtain that L is the unique minimizer of 



the functional 



ir{{A- LA)A-^{A- LAf)=tr{{I - L)A{I - L^)) , L e Sei , 



i.e., of the functional (3.17). From the proof of Lemma |3.2| , with B = A, it follows 
that the minimization problem 



min tr((/-L)A(/-L^)) 



decouples into seperate minimization problems (cf. ( |3.8| )) for the rows of L: 

min {-2lf + if Ali} 



(3.20) 
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for all i with > 0. Here ij and aj are the ith rows of L and A, respectively. The 
minimization problem corresponding to (3.18) is 



min W((I - L)A{I - L'^))u = min W(Au ~ 2lJ a, + if Ak) 



i=l 



This decouples into the same minimization problems as in (B.20). Hence the func- 
tionals in (3.17) and ( |3.18 ) have the same minimizer. 



Let J — diag((/ — L)A{I — L^)). Using the construction of G^i in (3.9a) we 
obtain 

DlJu = {D{,I - L)A{I - L^)D)u = (G^^^G^O^ 



'^^{G El A)ik(G Ei)ik — ^ Sik{GE' 



I )ik 



fc=l 



fe=l 
{i,k)eE' 



(G 



E' ) ii 



Dr. 



Hence Da — J^j ^ holds for all i, i.e., ( 3.19| ) holds. □ 

Corollary 3.6. From ( |3.19| ) it follows that diag(GE!^G£;!) = diag(G£;i) holds 
and thus, using ( 3.9b| ) we obtain 



diag(G£;z AG El 



(3.21) 



for the sparse approximate inverse Gei ■ 

The following theorem gives a main result in the theory of approximate inverses. 
It was first derived in [0|. A proof can be found in 0, too. 



Theorem 3.7. Let Gei be the approximate inverse in (S.i). Then Gei is the 
unique minimizer of the functional 



G 



T^trjGAG^ 
det (GAG'^)l 



G ^ Sei ■ 



(3.22) 



Proof. For G G Sei we use the decomposition G = D{I — L), with D diagonal 
and L G Se' ■ Furthermore, for L G 5^;! , Jl ■— diag((/ — L)A{I — L^)). Now note 



itr(GAG^) 
det (GAG^)^ 



det(G2)^ 



det(A)- 



det{D^JL) 



rdet(JL)" > det(A)-" det(Ji)' 



(3.23) 



The inequality in (3.23) follows from the inequality between the arithmetic and geo- 
metric mean: i J2i=i ^ (DLi "0^^" for > 0- 



For Gei in ( 3.9a ) we use the decomposition Ge' = D{I~L). For the approximate 



inverse Ge^ we then have Gei = (diag(G£;i)) "^Gei = D^i^I — L). From Lemma 
( |3.18 ) it follows that det( Jl) > det( J|_) for all L G Se' ■ Furthermore from Lemma 



( p. 19 
(E23 



3^ 
3.5 



we obtain that for Gei — D^{I ~L) we have [D^ Y^t ^ ^ ^"^^ t\ms equality in 
for G — Gei ■ We conclude that G^i is the unique minimizer of the functional 
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(HI)- 



Remark 3.8. The quantity 



K{A) 



dct(A)^ 



can be seen as a nonstandard condition number (cf . ^ ) . Properties of this quantity 
are given in H (Theorem 13.5). One elementary property is 



1 < K{A) < ^ = n{A) 



□ 



Corollary 3.9. For the approximate inverse G^i as in (3.9) we have (cf.(3.21)) 

\<K{GeiAGI,) ^ 



Aet{GEiAGl,Y 



I.e., 



d{A) < det{G%)-i = l[(GE');r = II^Ge^)-^ 



(3.24) 



i=l 



Let & be a lower triangular sparsity pattern that is larger than i?', i.e., E'- C & C 
{{hj) I 1 < j < * < n}. From the optimality result in Theorem [3.7] it follows that 

1 < KiGE:AG%) < KIGeiAGI,) . (3.25) 

□ Moti- 
vated by the theoretical results in Corollary 3.£ we propose to use the sparse approxi- 



mate inverse Ge' as in (3.9) for approximating d{A): Take d{GEi) = d{GEi) as 



an estimate for d{A). Some properties of this method are discussed in the following 
remark. 

Remark 3.10. We consider the method of approximating d{A) by d{GEi)^'^ = 
d{GE')~^- The practical realization of this method boils down to chosing a sparsity 
pattern and solving the (small) systems in ( 3.11 ). We list a few properties of this 
approach: 

1. The sparse approximate inverse exists for every symmetric positive definite 
A. Note that such an existence result does not hold for the incomplete Cholesky 
factorization. Furthermore, this factorization is obtained by solving low dimensional 



symmetric positive definite systems of the form PiAP^ gi = (cf. (3.11)), which can 
be realized in a stable way. 

2. The systems PiAP^ gi — Bi, i — 1,2, ... ,n, can be solved in parallel. 

3. For the computation of d{GE')^'^ — d{GE')~^ we only need the diagonal 
entries of Ge' (cf- (3.24)). In the systems PiAP^gi — Bi we then only have to 
compute the last entry of gi, i.e. {gi)ni ■ If these systems are solved using the Cholesky 
factorization, PiAP^ LiLf {Li lower triangular) we only need the {ni,ni) entry 
of Li, since (5i)„, = {Li)-f^^^. 

4. The sparse approximate inverse has an optimality property related to the 
determinant: The f uncti onal G K{GAG'^) , G G S e' , is minimal for Ge'- From 
this the inequality ([}.24|) and the monotonicity result (3.25) follow. 
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5. From ( ^^2^ ) we obtain the upper bound for the relative error d{A) / d{G ^i)^'^ — 
1. In §^we win derive useful lower bounds for this relative error. These are a posteriori 
error bounds which use the matrix Ge' ■ ^ 

4. A posteriori error estimation. In the previous section it has been ex- 
plained how an estimate d{GEi)^^ of d{A) can be computed. From ( ^.24| ) we have 
the error bound 



diA) 



d{GE^)- 



< 1 . 



(4.1) 



In this section we will discuss a posteriori estimators for the error d{A)/d{GEi) ^■ 
In §4.1 we apply the analysis from O] to derive an a posteriori lower bound for the 



quantity in (4.1). This approach results in safe, but often rather pessimistic bounds 
for the error. In §4.2 we propose a very simple stochastic method for error estimation. 



This method, although it does not yield guaranteed bounds for the error, turns out 
to be very useful in practice. 

4.1. Error estimation based on bounds from . In this section wc show 
how the analysis from ||^ can be used to obtain an error estimator. We first recall a 
main result from |^ (Theorem 2). Let ^ be a symmetric positive matrix of order n, 
m = tr(^), fj,2 = \\A\\l. and a{A) C [a, l3] with a > 0. Then: 



exp ( —[In a IntA 
\ n 



exp — [ln/3 In tu] 
\ n 



a 


tl ' 


-1 




a" 


''I 






a 




-1 



















< diA) < 



(4.2) 



where U - ^ 
In f[ 



this result is applied to obtain computable bounds for d{A). Often these 
bounds yield rather poor estimates of d(A). In the present paper we approximate 
d{A) by d(GE')~'^ and use the result in (4.2) for error estimation. The upper bound 
(4.1) turns out to be satisfactory in numerical experiments (cf. Therefore we 

restrict ourselves to the derivation of a lower bound for d{A) / d{G ei )~'^ , based on the 
left inequality in (4.2). 

Theorem 4.1. Let G^i be the approximate inverse from (3^) and < a 



(Ge'AG^i), ^ 

and 



IGe'AG 



E' 



< 



/i — 1 . The following holds: a < 1, 6 > 



exp 



1 



(51na + (l-a)2ln(l + 



1 



< 



diA) 



diGEi)- 



< 1 



(4.3) 



is already given in (4.1). We introduce the 
for the eigenvalues of Ge'AG^,. From (3.21) we obtain 



Proof. The right inequality in 
notation ti < r2 < . . . < r, 

i X)"=i '''i — ^ ^'^d from this it follows that a < ti < 1 holds. Furthermore, 
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yield s ^ — itr((G'^! AG^, )^) > 1 and thus (5 > 0. We now use the left inequality in 
(4.2) applied to the matrix G^iAG^, . Note that 



/Lti = tr(GBiAG|i) = n, fi2 n^, U 
A simple computation yields 



a^J■l - ^J■2 
an — 111 



1 - a 



a ti 


-1 


Ml 


1 


s 






. 


ti — a 


1 - a _ 



and 



ti - a 



1 - a 



(4.4) 



(4.5) 



Substitution of (4.5) in (4.4) results in 
— [In a Inti] 



a ti 


-1 


Ml 






. ^2 



1 



1 

(l-a)2+5 



(5\na+ (1 -aflnti) 
(5 In a + (1 - a f\n{l + 



1 - a 



Using this the left inequality in (^^) follows from the left inequality in (4.2). 



Note that for the lower bound in ( |4.3| ) to be computable, we need /i — ^\\G AG'^, 
and a strictly positive lower bound a for the smallest eigenvalue of G^i^G^, . We now 
discuss methods for computing a and fi. These methods are used in the numerical 
experiments in §^ 

We first discuss two methods for computing a. The first method, which can be 
applied if A is an M-matrix, is based on the following lemma, where we use the 
notation 1 = (1, 1, . . . , 1)^ e M". 

Lemma 4.2. Let A be a symmetric positive definite matrix of order n with Aij < 
for all i ^ j and G^i its sparse approximate inverse ). Furthermore, let z be such 
that 



WGpiAGl.z 



< 77 < 1 . 



Then 



holds. 



il-rj)\\z\\^' <X^UGE'AGlr, 



(4.6) 



Proof. From the assumptions it follows that A is an M-matrix. In g (Theorem 
4.1) it is proved that then G^'AG'^i is an Af-matrix, too. Let z* = (G^iAG^^i)'^!. 
Because {Ge'AG^i)~^ has only nonnegative entries it follows that 

\\{GeiAGI,)-^\\oo = \\z*\\oo = \\z+{z*-z)\\^ 

< ||z|U + \\{GEiAGl,)-^\U\GEiAGl,z ~ 1||^ 

< \\z\\^ + \\{GE^AGl,r'\\ooV ■ 
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Hence UG El AGl,)-^^ > {l~v)\\4^'- Using A„,i„( G^MGi;, ) > ||((G£,AG^,) 



we obtain the result (4.6). □ 

Based on this lemma we obtain the following method for computing a. Choose 
< ry ^ 1 and apply the conjugate gradient method to the system Ge'AG'^,z* = 1. 
This results in approximations z^, ... of z* . One iterates until the stopping crite- 
rion := \\GeiAGI:z' - Ijloo < ?7 is satisfied. Then take a (1 - d^)\\z^ . In 
view of efficiency one should not take a very small tolerance rj. In our experiments 
in ^ we use rj — 0.2 and = 1. Note that the CG method is applied to a system 
with the preconditioned matrix Gei AG'^i . In situations where the preconditioning is 
effective one may expect that relatively few CG iterations are needed to compute 
such that \\Ge'AG'^,z^ — 1||oo < ?7 is satisfied. Results of numerical experiments are 
presented in §||. 

As a second method for determining a, which is applicable to any symmetric positive 
definite A, we propose the Lanczos method for approximating eigenvalues applied to 
the matrix G^i AG^i . This method yields a decreasing sequence 
^ > ^miniG El AG^i) of approximations A^^ of Ainin(G£;i^G^;). If 

X[''> -\^in{GE'AGl,)<e (4.7) 



holds, then a ~ A^"*' ~e can be used in Theorem 4.1. However, in practice it is usually 



not known how to obtain reasonable values for e in (4.7). Therefore, in our experi- 



ments we use a simple heuristic for error estimation (instead of a rigorous bound as 



in (4.7)), based on the observed convergence behaviour of Ap'' (cf. §^). 
It is known that for the Lanczos method the convergence to extreme eigenvalues is 
relatively fast. Moreover, it often occurs that the small eigenvalues of the precon- 
ditioned matrix G^'^G^i are well-separated from the rest of the spectrum, which 

has a positive effect on the convergence speed A^''^ — > XminiG e' AG^i) . In numerical 
experiments we indeed observe that often already after a few Lanczos iterations we 
have an approximation of Xi-ain{G e^ AG^i) with an estimated relative error of a few 
percent. However, for the a computed in this second method we do not have a rig- 
orous analysis which guarantees that < a < Xnun{G e' AG'^,) holds. Nevertheless, 
from numerical experiments we see that this method is satisfactory. This is partly ex- 
plained by the relatively fast convergence of the Lanczos method towards the smallest 



eigenvalue. A further explanation follows from the form of the lower bound in (4.2) 



For a <S 1, S ^ 1, which is typically the case in our experiments in i|q, this lower 
bound essentially behaves like exp(i51nQ;) =: 5(a). Note that < ^-^^y- — S <^ I 
holds. Hence the sensitivity of the lower bound with respect to perturbations in a is 
very mild. 

We now discuss the computation of the quantity /i — ^\\Ge'AG'^i \ \e, which is needed 



in (4.3). Clearly, for computing /i one needs the matrices G^i and A. To avoid un- 
necessary storage requirements one should not compute the matrix X := Ge'AG'^i 
and then determine /x = i||Ar|||,. A with respect to storage more efficient approach 
can be based on 

n 

\\GEiAGUl^Y.\\^E'AGl,e,\\l , 

i=l 

where is the ith basis vector in M". For the computation of ||G£;i^GUiei||2, 
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i = 1,2, .. . , n, which can be done in paraUel, one needs only sparse matrix- vector 
multiplications with the matrices G^i and A. Furthermore, for the computation 
of AG^iCi one can use that {DGEiA)ij — {G E'^)ii — ^ij for G E'', with 

D := diag{G El) (cf. i^). It follows from (|3?10| ) that 



AGl.e, = {I- P^Pi)AGl,e, + P,AGl^e, 

= {I- PfP,)AGlie, + PfP,AGl,D-^ei 
^{I-P^Pi)AGl,e,+D^^e, 

holds. 

Remark 4.3. Note that for the error estimators discussed in this section the 
matrix G^i must be available (and thus stored), whereas for the computation of 
the approximation d{GEi)~^ of d{A) we do not have to store the matrix Gei (cf. 
Remark 3.1C item 3). Furthermore, as we will see in §|], the computation of these 



error estimators is relatively expensive. □ 

4.2. Error estimation based on a Monte Carlo approach. In this section 
we discuss a simple error estimation method which turns out to be useful in practice. 
Opposite to those treated in the previous section this method does not yield (an 
approximation of) bounds for the error. 
The exact error is given by 



diGE'Y 



diGEiAG'Ei)=di£E') , 



where £ei ■— GeiAG^i is a sparse symmetric positive definite matrix. Fore ease of 
presentation we assume that the pattern is sufficiently large such that 

p{I~£ei)<1 (4.8) 



holds. In |ll[] it is proved that if A is an M-matrix or a (block) i/-matrix then (4.8) 
is satisfied for every lower triangular pattern E^ . In the numerical exp eriments (cf. 
§1^) with matrices which are not Af -matrices or (block) j J-ma trices ( |4.8[ ) turns out to 
be satisfied for standard choices of E^ . We note that if ( [4.8[) does not hold then the 
technique discussed below can still be applied if one replaces £ei by ijj£ei with w > 
a suitable damping factor such that p{I — uj£ei ) < 1 is satisfied. 

For the exact error we obtain, using a Taylor expansion of ln(/ — B) for B G 
with p{B) < 1 (cf. |]): 

d{£Ei) ~ exp ( —\n{det{£ El))] = exp ( —tr{\n{£Ei)) 
\n I \n 



pn X n 



exp (^^tr(ln(/ - (/ - £^0))) - exp (-^t^E f^'^ )j (4-9) 
Hence, an error estimation can be based on estimates for the partial sums Sm := 



k 



X^fcLi itr((/ — f^;;)*^). The construction of Gei is such that diag(£'£i) = / (cf. (3.21)) 
and thus tr(f£i) = n and Si — 0. For 5*2 we have 



itr((/ - £ei )') - itr(/ - 2£ei +£l)^-^n+^ 



(4.10) 
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For 5*3 we obtain 

D / .3 (4.11) 

Note that in 5*2 and 5*3 the quantity tr(g^;) = H^biIIf — ll^^i AG^i |||. occurs which 
is also used in the error estimator in |4.l| . In this section we use a Monte Carlo 
method to approximate the trace quantities in Sm- The method we use is based on 
the following proposition ||, ||] . 

Proposition 4.4. Let H be a symmetric matrix of order n with tr(iJ) ^ 0. 
Let V be the discrete random variable which takes the values 1 and —1 each with 
probability 0.5 and let z be a vector of n independent samples from V . Then Liz is 
an unbiased estimator oftr{H): 

E{z'^Hz) = tr{H) , 

and 

var{z'^Hz) = 2 ^ /i^ . 

For approximating the trace quantity in iS'2 we use the following Monte Carlo algo- 
rithm: 



for j = l,2,...,M 

1. Generate Zj € M" with entries which are uniformly distributed in (0, 1). 

2. If {zj)i < 0.5 then {zj)i —1, otherwise, {zj)i 1. 

3. yj := Se'Zj, ^j-.^yjyj. 



Based on Proposition 4.4 and (4.10) we use 

1 



^2 



—n 

2 2M 



1 



(4.12) 



as an approximation for ^2. The corresponding error estimator is 

E2 = exp(--S'2). 
n 

For the approximation of S'3 we replace step 3 in the algorithm above by 

i.yj-.^EEiZj, yj-.^EEiyj, -.^ ^yj yj - \yj yj 
and we use 



(4.13) 



7 1 ^ 



(4.14) 



as an estimate for S'3. The corresponding error estimator is 

n 



(4.15) 
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Clearly, this technique can be extended to the partial sums Sm with m > 3. However, 
in our applications we only use S2 and for error estimation. It turns out that, at 



least in our experiments, the two leading terms in the expansion (4.9) are sufficient 
for a reasonable error estimation. Note that due to the truncation of the Taylor ex- 
pansion, the estimators E2 and Eg, for d{£^i ) are biased. 

It is shown in ^ that based on the so-called Hoeffding inequality (cf. [|l3|) prob- 
abilistic bounds for |-p J2iLi zj Hzi — tr(_ff)| can be derived, where zi, Z2, ■ ■ ■ , zm are 



independent random variables as in Proposition iA . In this paper we do not use 
these bounds. Based on numerical experiments we take a fixed small value for the 
parameter M in the Monte Carlo algorithm above (in the experiments in §5: M — Q). 



Remark 4.5. In the setting of this paper Proposition 4.4 is applied with H = 



p{£^i), where p is a known polynomial of degree 2 or 3. In the Monte Carlo technique 



for approximating det(yl) — exp(tr(ln(A))) from Proposition 4.4 is applied with 
H = ln{A). The quantity z'^\n{A)z, which can be considered as a Ricmann-Stieltjes 
integral, is approximated using suitable quadrature rules. In |^ this quadrature is 
based on a Gauss-Christoffcl technique where the unknown nodes and weights in the 
quadrature rule are determined using the Lanczos method. For a detailed explanation 
of this method we refer to [|. 

A further alternative that could be considered for error estimation is the use of this 
method from j^. In the setting here, this method could be used to compute a (rough) 
approximation of det(G£;iylG'^,)^/". We did not investigate this possibility. The re- 
sults in ^ give an indication that this alternative is probably much more expensive 
than the method presented in this section. □ 

5. Numerical experiments. In this section we present some results of numer- 
ical experiments with the methods introduced in ^ and All experiments are done 
using a MATLAB implementation. We use the MATLAB notation nnz{B) for the 
number of nonzero entries in a matrix B. 



Experiment 1 (discrete 2D Laplacian). We consider the standard 5~-point discrete 
Laplacian on a uniform square grid with m mesh points in both directions, i.e. n = m? . 
For this symmetric positive definite matrix the eigenvalues are known: 

A.^ = 4(to + If Uv?{ ) + sin^( ^^^j ) , l<v,^i<m. 

\ 2[m + Ij 2(m + 1) J 

For the choice of the sparsity pattern we use a simple approach based on the 
nonzero structure of (powers of) the matrix A: 

EHk):^{{i,j)\i>j andiA%^0} , fc = l,2,.... (5.1) 

We first describe some features of the methods for the case m = 30, k — 2 and after 
that we will vary m and k. Let A denote the discrete Laplacian for the case m = 30 and 
La its lower triangular part. We then have nnz{LA) = 2640. For the sparse approxi- 
mate inverse we obtain nnz{G {2)) — 6002. Th e sys tems PiAP^gi — (0, 0, . . . ,0, 1)"^ 
that have to be solved to determine Gei(2) (cf. ( ^.ll| )) have dimensions between 1 and 
7; the mean of these dimensions is 6.7. As an approximation of d{A) = 3.1379 10'^ we 
obtain 

n 

rf(Gi5<(2))"' = d{GE^(2))-' = n(GiJ'(2))i^ " = 3.2526 10^ . 

1=1 
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Hence d( A)/ (1(0^1 [2)) ^ = 0.965. For t he computation of this approximation along 
the hnes as described in Remark 3.10| , item 3, we have to compute the Cholesky 



factorizations PiAP^ = LiLf , i = 1,2, .. . ,n. For this approximately 41 10'^ flops 
are needed (in the MATLAB implementation) . If we compare this with the costs of 
one matrix-vector multiplication A* x (8760 flops), denoted by MATVEC, it follows 
that for computing this approximation of d{A), with an error of 3.5 percent, we need 
arithmetic work comparable to only 5 MATVEC. 

We will see that the arithmetic costs for error estimation are significantly higher. 
We first consider the methods of |4.l| . The arithmetic costs are measured in terms 
of MATVEC. For the computation of a as indicated in Lemma |4.2| with 77 — 0.2, 
using the CG method with starting vector 1 = (1, 1, ... , 1)^ we need 8 iterations. 
In each CG iteration we have to compute a matrix-vector multiplication G^iAG^j^iX, 
which costs approximately 3.7 MATVEC. We obtain acG = 0.0155. For the method 
based on the Lanczos method for approximating Xj^iniG e' ^G^i) we use the heuristic 
stopping criterion 

\X[''> -X['-'^\<Qm\x[''>\ . (5.2) 

We then need 7 Lanczos iterations, resulting in ckLanczos = 0.0254. A direct computa- 
tion results in Xn^iniG e' AG'^^) = 0.025347. 

For the computation of /i = WG^iAG'^i |||, we first computed the lower triangular part 
oi X — G^iAG'^i and then computed ||A^||_f (making use of symmetry). The total 



costs of this are approximately 18 MATVEC. Application of Lemma 4.1, with acG 
and aLanczos yields the two intervals 

[0.880, 1] and [0.895, 1] , 

which both contain the exact error 0.965. In both cases, the total costs for error 
estimation are 40-45 MATVEC, which is approximately 10 times more than the costs 
for computing the approximation d{G ^u^^)''^ ■ 

We now consider the method of §[4 2[ We use the estimators E2 and from 



( [4. 131) , ( [4.15D with M = 6. The results are E2 = 0.980, Ey. = 0.973. Note that 
the order of magnitude of the exact error (3.5 percent) is approximated well by both 
E2 (2.0 percent) and E^ (2.7 percent). In step 3 in the Monte Carlo algorithm 
for computing >S'2 we need one matrix- vector multiplication G^iAG^'^iX (costs 3.7 
MATVEC). The total arithmetic costs for E2 are approximately 20 MATVEC. For 
we need two matrix-vector multiplications with S^i in the third step of the Monte 
Carlo algorithm. The total costs for E^ are approximately 40 MATVEC. 

In Table |5T| we give results for the discrete 2D Laplacian with m = 30 (n = 900), 
TO = 100 (n = 10000) and m = 200 (n = 40000). We use the sparsity pattern E\2). 
In the third column of this table we give the computed approximation of d{A) and the 
corresponding relative error. In the fourth column we give the total arithmetic costs 



for the Cholesky factorization of the matrices PiAP^ , i = 1,2, . . . ,n (cf. Remark 3.1C , 
item 3). In the columns 5-8 we give the results and corresponding arithmetic costs 
for the error estimators discussed in §0. The fifth column corresponds to the method 



discussed in §4.1 with a determined using the CG method applied to G^iAG^i 



with starting vector 1. In the stopping criterion we take 77 = 0. 2 (c f. Lemma 4.2) 



The computed a = acG is used as input for the lower bound in (L3). The resulting 
bound for the relative error and the arithmetic costs for computing this error bound 
are shown in column 5. In column 6 one finds the computed error bounds if a is 



determined using the Lanczos method with stopping criterion (5.2). In the last two 
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Table 5.1 

Results for 2D discrete Laplacian with = (2) 



n. 


d{A) 


'^('^B'(2)) ^ 


costs for 


Thm. [1.1|, 


Thm. 


MC 


MC 






(error) 


'^(G'e'(2))~^ 


"CG 


Q^Lanczos 


E2 


Ez 


900 


3.138 10^ 


3.253 10^ 


5 MV 


< 12% 


< 11% 


2.0% 


2.7% 






(3.5%) 




(45 MV) 


(45 MV) 


(20 MV) 


(40 MV) 


10000 


3.292 10" 


3.434 10'' 


5 MV 


< 21% 


< 19% 


2.2% 


2.6% 






(4.1%) 




(140 MV) 


(102 MV) 


(24 MV) 


(48 MV) 


40000 


1.300 10^ 


1.359 10^ 


5 MV 


< 26% 


< 24% 


2.2% 


2.7% 






(4.3%) 




(276 MV) 


(159 MV) 


(24 MV) 


(48 MV) 



Table 5.2 

Results for 2D discrete Laplacian with E^ = E^ (4) 



n 


d{A) 


(error) 


costs for 
'^(G'b'{4))~^ 


Thm. |t.l|, 


Thm. 

'^Lanczos 


MC 

E2 


MC 

E3 


900 


3.138 10^ 


3.177 10^ 
(1.2%) 


41 MV 


< 3.5% 
(137 MV) 


< 3.0% 
(146 MV) 


0.65% 
(54 MV) 


1.1% 
(108 MV) 


10000 


3.292 10* 


3.347 10* 
(1.6%) 


45 MV 


< 7.7% 
(263 MV) 


< 7.0% 
(226 MV) 


0.91% 
(55 MV) 


1.1% 
(110 MV) 


40000 


1.300 10^ 


1.323 10^ 
(1.7%) 


46 MV 


< 10% 
(487 MV) 


< 9.3% 
(348 MV) 


0.93% 
(56 MV) 


1.1% 
(112 MV) 



columns the results for the Monte Carlo estimators E2 (4:-13) and £'3 (4.15) are given. 
In Table 5.2 we show the results and corresponding arithmetic costs for the method 
with sparsity pattern — i?'(4). 

Concerning the numerical results we note the following. From the third and fourth 
column in Table 5T we see that using this method we can obtain an approxima- 
tion of d{A) with relative error only a few percent and arithmetic costs only a few 
MATVEC. Moreover, this efhciency hardly depe nds on t he d imension n. Comparison 
of the third and fourth columns of the Tables 5.1 



and 5.: 



shows that the approx- 
imation significantly improves if we enlarge the pattern from E^{2) to _E'(4). The 
corresponding arithmetic costs increase by a factor of about 9. This is caused by 
the fact that the mean of the dimensions of the systems PiAP^ , i = 1, 2, n, in- 
creases from approximately 7 {E^{2)) to approximately 20. For n — 10000 wc have 
nnz{LA) = 29800, nnz{GEi{2)) = 69002, nnz{GEi(i)) = 204030. For the other n 
values we have similar ratios between the number of nonzeros in the matrices La and 
G^i . Note that the matrix G^i has to be stored for the error estimation but not 
for the computation of the approximation (1(0^1)^"^ ■ The error bounds in the fifth 
and sixth column in the Tables 



5.1 



and p.2| are rather conservative and expensive. 
Furthermore there is some deterioration in the quality and a quite strong increase in 
the costs if the dimension n grows. The strong increase in the costs is mainly due to 
the fact that the CG and Lanczos method both need significantly more iterations if n 
increases. This is a well-known phenomenom (the matrix G^iAG^^i has a condition 
number that is proportional to n) . Also note that the costs for these error estimators 
are (very) high compared to the costs of the computation of d{GE')^^- The results 
in the last two columns indicate that the Monte Carlo error estimators, although less 
reliable, are more favourable. 

In Figure jA we show the eigenvalues of the matrix G^iAG'^i for the case n — 900, 
E^ — E\2) (computed with the MATLAB ftmction eig). The eigenvalues are in the 
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interval [0.025,1.4]. The mean of these eigenvalues is 1 ( 11(6^1 AG'^,) — 1). One 
can see that relatively many eigenvalues are close to 1 and only a few eigenvalues are 
close to zero. 



Fig. 5.1. Eigenvalues of the matrix G^iAG'^i in Experiment 1 




100 200 300 400 500 600 700 800 900 



Experiment 2 (MATLAB random sparse matrix). The sparsity structure of the 
matrices considered in Experiment 1 is very regular. In this experiment we con- 
sider matrices with a pattern of nonzero entries that is very irregular. We used the 
MATLAB generator (sPRAND(n, n, 2/n)) to generate a matrix B of order n with ap- 
proximately 2n nonzero entries. These are uniformly distributed random entries in 
(0, 1). The matrix B is then sparse symmetric positive semidefinite. In the generic 
case this matrix has many eigenvalues zero. To obtain a positive definite matrix we 
generated a random vector d with all entries chosen from a uniform distribution on 
the interval (0, 1) {d :=RAND(n, 1)). As a testmatrix we used A := i?-^_B + diag(d). We 
performed numerical experiments similar to those in Experiment 1 above. We only 
consider the case with sparsity pattern — E^{2). The error estimator based on the 
CG method is not applicable because the sign condition in Lemma 4.2 is not fulfilled. 
For the case n = 900 the eigenvalues of A and of G^'AG'^, are shown in Figure 5.2. 
For A the smallest and largest eigenvalues are 0.0099 and 5.70, respectively. The 
picture on the right in Figure shows that for this matrix A sparse approximate 
inverse preconditioning results in a very well-conditioned matrix. Related to this, one 
can see in Table 5.3 that for this random matrix A the approximation of d{A) based 
on the sparse approximate inverse is much better than for the discrete Laplacian in 
Experiment 1. For n = 900,10000,40000 we obtain nnz{LA) = 2730,29789,120216 
and nnz{GEi) = 7477,82290,335139, respectively. For n ^ 900,10000,40000 the 
mean of the dimensions of the systems PiAPf , i = 1,2,... , n is 10.6, 10.8, 11.0, re- 
spectively. In all three cases the costs for a matrix- vector multiplication GeiAGeix 
are approximately 4.3 MV. Furthermore, in all three cases the matrix Gei^AG^i is 
well-cond ition ed and the number of Lanczos iterations needed to satisfy the stopping 
criterion (5_^) hardly depends on n. Due to th is, f or increasing n, the growth in the 
costs for the error estimator base d o n Th eore m |4.l| (co lum n 5) is much slower than in 
Experiment 1. As in the Tables 5T and |]^, in Table 5^ the error quantities in the 
columns 3, 5,6,7 are bounds or estimates for the relative error 1 — d{GEiAG'^i). 
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Fig. 5.2. Eigenvalues of the matrices A and G^iAG^i in Experiment 2 




Table 5.3 

Results for MATLAB random sparse matrices with E^ = E^{2) 



n 


d{A) 


d{G,s^)-' 
(error) 


costs for 


Thm. 4.1, 

<^Lanczos 


MC 

E2 


MC 


900 


0.82453 


0.82521 
(8.3 10-'') 


23 MV 


< 9.8 10-"' 
(110 MV) 


1.4 10-3 
(26 MV) 


1.0 1.0-3 
(52 MV) 


10000 




0.81053 

(-) 


18 MV 


< 1.1 10-3 
(139 MV) 


8.4 10-4 
(26 MV) 


7.4 10-4 
(52 MV) 


40000 




0.82033 

(-) 


18 MV 


< 1.0 10-3 
(146 MV) 


6.2 lO-"* 
(26 MV) 


8.3 10-4 
(52 MV) 



For n = 10000,40000 the values of d{A) are not given (column 2). This has to do 
with the fact that for these matrices with very irregular sparsity patterns the Cholesky 
factorization A — LL^ suffers from much more fill-in than for the matrices in the Ex- 
periments 1 and 3. For the matrix A in this experiment with n = 900 we have 
nnz{LA) = 2730 and nnz(L) = 72766. For n = 10000 we run into storage problems 
if we try to compute the Cholesky factorization using the MATLAB function CHOL. 

Experiment 3 (QCD type matrix). In this experiment we consider a complex 
Hermitean positive definite matrix with sparsity structure as in Experiment 1. This 
matrix is motivated by applications from the QCD field. In QCD simulations the 
determinant of the so-called Wilson fermion matrix is of interest. These matrices 
and some of their properties are discussed in |^ ^ . The nonzero entries in a Wilson 
fermion matrix arc induced by a nearest neighbour coupling in a regular 4-dimensional 
grid. These couplings consist of 12 x 12 complex matrices M^y, which have a tensor 
product structure M^y = P^y ® Uxy, where P^y G R'*^'* is a projector, Uxy G C'^^"^ is 
from SU3 and x and y denote nearest neighbours in the grid. These coupling matrices 
Mxy strongly fluctuate as a function of x and y. Here we consider a (toy) problem 
with a matrix which has some similarities with these Wilson fermion matrices. We 
start with a 2-dimcnsional regular grid as in Experiment 1 (n grid points). For the 
couplings with nearest neighbours we use complex numbers with length 1. These 
numbers are chosen as follows. The couplings with south and west neighbours at a 
grid point x are exp(2«7ras(a;)) and exp(2i7ravK(a;)), respectively, where as{x) and 
awix) are chosen from a uniform distribution on the interval (0, 1). The couplings 
with the north and east neighbours are taken such that the matrix is hermitean. To 
make the comparison with Experiment 1 easier the matrix is scaled by the factor n. 
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i.e. the couplings with nearest neighbours have length n. For the diagonal we take 7/, 
where 7 is chosen such that the smallest eigenvalue of the resulting matrix is approx- 
imately 1 (this can be realized by using the MATLAB function EIGS for estimating 
the smallest eigenvalue). We performed numerical experiments as in Experiment 1 
with ~ E^{2). The number of nonzero entries in La and Gei are the same as 
in Experiment 1. F or n = 900 the eigenvalues of the matrices A and GeiAG^i are 
shown in Figure IsTs] . These spectra are in the intervals [1, 6.6 10^] and [1.7 10"^, 1.5], 
respectively. 



The results of numerical experiments are presented in Table 5.4. Note that the error 



Fig. 5.3. Eigenvalues of the matrices A and G^iAG^i in Experiment 3 





estimator from §4.1 in which the CG method is used for computing a can not be used 



for this matrix (assumptions in Lemma 4.2 are not satisfied). We did not consider the 
case n — 40000 here because then the application of the EIG function for computing 
the smallest eigenvalue led to memory problems. 



Comparison of the results in Table 5.4 with those in Table 5.1 shows that when the 



Table 5.4 
Results for QCD type matrix with 



E'{2) 



n 


d{A) 


(error) 


costs for 
d(GsO-' 


Thm. 4.1, 

"^Lanczos 


MC 
E2 


MC 
E3 


900 


2.500 10=^ 


2.620 10^ 
(4.6%) 


5 MV 


< 24% 
(79 MV) 


2.6% 
(23 MV) 


3.3% 
(46 MV) 


10000 


2.739 10"' 


2.842 10* 
(3.6%) 


5 MV 


< 31% 
(133 MV) 


2.4% 
(23 MV) 


2.7% 
(46 MV) 


22500 


6.173 10'' 


6.391 10* 
(3.4%) 


5 MV 


< 32% 
(198 MV) 


2.3% 

(23 MV) 


2.8% 
(46 MV) 



method is applied to the QCD type of problem instead of the discrete Laplacian the 
performance of the method does not change very much. 

Finally, we note that in all measurements of the arithmetic costs we did not take 
into account the costs of determining the sparsity pattern E^{k) and of building the 
matrices Pi APT ■ 
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